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Dissipation from harmonic energy eigenstates is used to characterize energy transport in binary 
isotopically disordered (BID) Fermi-Pasta-Ulam (FPU-/3) chains. Using a continuum analog for 
the corresponding harmonic portion of the Hamiltonian, the time-independent wave amplitude is 
calculated for a plane wave having wavelength A that is incident upon the disordered section, and 
the solution is mapped onto the discrete chain. Due to Anderson localization, energy is initially 
localized near the incident end of the chain, and in the absence of anharmonicity the wave amplitude 
is stationary in time. For sufficient anharmonicity, however, mode transitions lead to dissipation. 
Energy transport along the chain is quantified using both the second moment M of the site energy, 
and the number of masses contributing to transport, which was estimated from the localization 
parameter. Over the time scales studied, M increased linearly in time, yielding an effective transport 
coefficient G. At low and intermediate impurity concentration c, G(c) can be characterized by a 
competition between the rate of mode transitions and the time for energy to propagate a distance 
equal to the localization length £. At the highest concentrations (1.6 < cA < 16.0), there is 
significant mode transition suppression in BID systems, and the transport coefficient G(c) becomes 
proportional to £(c). 

PACS numbers: 45.05.+X, 45.30.+S, 46.40.Cd, 62.20.Pw, 



I. INTRODUCTION 

Nonlinear binary disordered chains are a useful sys- 
tems for studying the essential characteristics of energy 
dissipation and transport in materials. Although binary 
disorder is an idealized model, it has practical application 
to a number of fields: isotopic disorder effects line width 
broadening in spectroscopyfi*2^ the glass transition has 
been considered in terms of binary changes in elasticity 
parameters^ and isolated mechanical defects that occur 
in practice can lead to a variety of nonlinear effectsi^ 

An interesting behavior of binary isotopic disorder 
(BID) occurs in discrete lattices, where the system un- 
dergoes a pure-disordered-pure transition as the impu- 
rity concentration varies from zero to one. For harmonic 
one-dimensional chains composed of discrete elements, fi- 
nite disorder destroys spatial invariance and gives rise to 
Anderson^ localization, characterized by spatially local- 
ized eigenstates. Given a discrete system in a localized 
energy eigenstate, the addition of anharmonicity will lead 
to interactions that create new modes that can propa- 
gate through the system. These propagating modes ei- 
ther will become localized or will undergo further mode 
transitions. 

We are interested in the rate of energy dissipation 
from localized disturbances. Previous studies of dis- 
sipation have used either a singular (one or very few 
elements) pulse or a Gaussian pulse for the initial 
displacement i^iiSiiiii^ Because neither is an eigenstate 
of the system, the singular pulse will spread and the 
Gaussian pulse will will propagate ballistically in both 



anharmonic and harmonic systems. For Gaussian and 
'kicked' singular pulses, the initial behavior is ballistic. 
In time, scattering and localization slow the rates of prop- 
agation and dispersion. 

To eliminate the initial ballistic motion, the system can 
be started in an energy eigenstate of the correpsonding 
harmonic chain. The middle section of the chain con- 
tains disorder, and there are "pure" sections at both ends. 
Given the location of each impurity in the disordered sec- 
tion of the chain, a solution is found for the continuum 
analog, with the boundary condition of an incident plane 
wave with frequency lo; there is a reflective wave and a 
transmitted transmitted wave. Because the continuum 
calculation is performed for the harmonic system, the so- 
lution is separable into spatial and temporal components. 
The time dependence is sinusoidal, and the entire system 
has a constant temporal phase. The continuum solution 
for the wave amplitude everywhere is mapped onto the 
discrete chain, and becomes the initial displacement. 

The advantage of the this approach is that in the ab- 
sence of anharmonicity the wave amplitude is stationary. 
For the harmonic chain, the wave remains localized indef- 
initely, and there is no energy transport. The addition of 
anharmonicity will give rise to mode transitions that will 
de-localize the wave and lead to energy transport through 
the chain. Instead of having a ballistic-diffusive transi- 
tion, this initial condition leads to a localized-diffusive 
transition. 

The remaining question is whether the energy will dis- 
sipate in a diffusive manner. In one-dimensional systems, 
mode interactions can lead to frequencies arbitrarily close 
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to zero. Because the localization length of these modes 
is proportional to uj~ 2 and the scattering cross section 
is proportional to w 2 , these low frequency modes should 
propagate ballistically in a finite system. Because the en- 
ergy content of these modes is proportional to uj 2 , how- 
ever, it is unlikely that these ballistic modes contribute 
substantially to energy transport. 

To study this energy transport, the following numerical 
experiment will use the binary disordered Fermi-Pasta- 
UlamH chain with quartic spring potentials (FPU-/?). 
Starting from the initial energy eigenstate, numerical in- 
tegration will be used to calculate the spatial distribu- 
tion of energy as a function of time. Based on a local 
concept of thermal transport^ an effective transport co- 
efficient will be sought from the second moment of the 
site energy, and the method is compared to the Helfand 14 
moments for thermal conductivity. The second moment 
will exhibit diffusive behavior, a fact that will be cor- 
roborated qualitatively from the number of masses that 
energy is distributed, which is estimated from the local- 
ization parameter 

II. NUMERICAL EXPERIMENT 
A. FPU-/3 Chain 



elsewhere |£»i£ was chosen as a compromise between speed 
and accuracy, and the results were insensitive to two-fold 
changes in At. Using this time step, the energy fluctua- 
tions were always less then 0.2 % (See Appendix). 



B. Semi-Infinite Approximation 

In most cases, low frequency waves will propagate 
down the chain and reach the far end. If the wave is 
reflected, it could affect the accuracy of the transport 
calculation. To mitigate this effect, 10 % of the masses 
at the far end were given a viscous force F V i S : 

F vis =-r)q (3) 

This approach has been used elsewhere to achieve a sim- 
ilar effect^ For these calculations, a viscosity r\ of 0.2 
was sufficient to eliminate the effects of reflection. 

The viscous damping, combined with the accurate time 
integration, simplified the task of identifying finite size 
effects. As the system length decreased and the local- 
ization length increased, the likelihood of a considerable 
amount of energy reaching the far end of the system in- 
creased. Fortunately, this occurrence was easily identified 
by changes in the system total energy. 



The FPU-/3 chain is composed of masses interacting 
with nearest neighbors through springs. The Hamilto- 
nian H of a chain having N masses is function of the 
mass momenta p, and the mass displacements about 
their equilibrium position: 
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The harmonic spring force coefficient K is set equal to 
one, and in the absence of impurities, each mass has the 
same value m =l. The site energy Ei is the sum of the 
kinetic energy plus one-half of the neighboring spring po- 
tential energies. The total energy Et is the sum of the 
site energies. 

Whenever possible, the results are expressed in dimen- 
sionless units through the use of appropriate scaling fac- 
tors. Time is scaled by the natural frequency lu of a 
single harmonic oscillator: 
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Lengths are scaled by the equilibrium mass separation 
distance a. 

The time-dependent behavior was determined by nu- 
merical integration using the sixth-order Yoshida 17 sym- 
plectic integration algorithm. Specifically, best results 
were obtained from the "Solution A" coefficients (see Ta- 
ble 1 in Ref . ITtT) . For the systems studied here, the time 
step At was approximately 1/200 the period of the ini- 
tial mode. This time step was consistent with that used 



C. Initial Displacement 

The initial condition was a stationary state exhibiting 
Anderson localization for the harmonic component of the 
Hamiltonian. This initial condition was chosen so that for 
(3=0 there is no net energy transport, and these systems 
could be used as a test to confirm the accuracy of the 
model and the numerical integration. 



The continuum Kronig-Penney liquid mode 



20.21 



IS 

used to estimate the initial displacement. For the har- 
monic FPU chain, the analogous continuum system is an 
elastic rod having mass density fx = m /a and Youngs 
modulus Y = Ka. Between the impurities, a longitudi- 
nal wave ip(x,t\u) with frequency u> will propagate with 
velocity q = \fYj Jl. A harmonic oscillator impurity, ap- 
proximated by a point defect, located at x 1 will give rise 
to a reactive force due to the impurity impedance Z to a 
wave with frequency lu: 22 - 
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For a harmonic system such as this, one can assume a 
sinusoidal solution with frequency u> (ip = <j>(x\u>) e~ luJt ): 
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where £ = oj/q. 
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Equation JSJ is solved for a system having impurities 
at integer locations with probability c. The boundary 
conditions are a unit amplitude incident wave and a re- 
flected wave at one end and only a transmitted wave at 
the other end. The initial displacement for the physi- 
cal system of masses and springs is taken from the real 
component of the solution (j){x\uj). 



D. Impurity Cross Section 

To put some of the results in a familiar context, it will 
be useful to characterize an impurity by its cross section 
to the original wave. The cross section tr of an individual 
impurity can be expressed as a function of the impurity 
impedance Z^ 



\Z(uj)\ 2 +4Km 

The isotopic impurities consist of a constant mass m+ 
added to the existing mass m ; m + may be either positive 
or negative. (The impurity mass m/ = m a + m + .) The 
impedance Z of a mass impurity in a continuum system^ 
is modified by c s , which is the ratio of the group velocity 
(v g = dui/dk) to the longitudinal velocity 

Z(u>) — —iujm + /c s (7) 

For the discrete chain, c s — cos(fca/2), where k is the 
wavenumber (2n/X) and A is the displacement wave- 
length. 



E. Continuum-Discrete Mapping 

The boundary condition for the chain is zero dis- 
placement at each end. Because the continuum solution 
(j)(x — 0, L\lu) will (with almost certainty) not equal zero 
at (x = 0, L), a method is needed for adjusting the con- 
tinuum solution to accommodate the constraints of the 
discrete chain. The continuum solution is first mapped to 
the discrete chain with no modification, and then, start- 
ing at the end and searching along the chain, the end 
is relocated at the mass having the smallest oscillation 
amplitude. The end is relocated to this mass, and its 
displacement is fixed at zero. The process is repeated at 
the opposite end of the chain. 

To facilitate this task, a suitable initial displacement 
wavelength A is needed to ensure that some mass has 
an equilibrium displacement acceptably close to zero. If 
the displacement wavelength is an integer multiple of the 
equilibrium mass displacement a, the oscillation ampli- 
tude of masses along each successive wavelength remains 
constant, until the next impurity changes the phase. 

For this experiment, a nominal displacement (equal to 
an integer multiple of the unit spacing) is chosen first. 
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FIG. 1: Initial displacements qi for a particular system having 
wavlength 31.8, impurity concentration 0.01, and impurity 
scattering cross section 0.5. Mass displacement are denoted 
by small circles and the impurity locations are denoted by 
large circles. The first impurity is located at i = L, and the 
dashed line is proportional to e - ' 1- '^. 



The working wavelength is the nominal wavelength mi- 
nus 0.2 unit spacings. With this wavelength, the dis- 
placement closest to zero repeats every 5 wavelengths, 
with 9 nodal points occurring at different fractions of the 
spacing a. Probabilistically, this reduces the minimum 
mass oscillation amplitude by nearly an order of magni- 
tude over an integer wavelength. 

An example input displacement is shown in Fig. ^ for 
a system with length 10 000, displacement wavelength 
31.8, impurity cross section 0.5, and impurity concentra- 
tion 0.01. The small dots in the figure represent the dis- 
placement of the masses. The larger filled circles denote 
the location of impurities. The effect of the impurities 
is to change both the the displacement amplitude and 
phase. Also shown in the figure is a dashed line that is 
proportional to e~^~ L ^^, where l is the location of the 
first impurity. Although this particular initial condition 
would suggest that the displacement amplitudes decrease 
monotonically, that is not always the case. 



F. Dense Systems 

For BID systems, the localization length is not a mono- 
tonic function of impurity concentration. Rather, the 
localization length has a minimum with respect to impu- 
rity concentration; at higher concentrations, the system 
returns to a "pure" system. The localization length £ 
of BID systems at dilute impurity concentrations, such 
that cA <C 1, can be calculated from the resistivity scal- 
ing lawt^ 

Ci = cln(l + p) (8) 
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FIG. 2: Initial displacements for a particular system hav- 
ing wavlength 31.8, impurity concentration 0.50, and all three 
impurity scattering cross sections. Mass displacement are de- 
noted by line, and the impurity locations are denoted by cir- 
cles. The systems were shifted horizontally so that the first 
impurity is located at i =200. 



The single impurity resistitivity p is related to the single 



impurity cross section a 
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As the impurity concentration increases, the system 
approaches a homogeneous system. In the limit c — ► 1, 
the system is again ordered, and the localization length 
diverges: 



Cii = (l-c)ln(l+p') 



(10) 



The quantity p' characterizes a system having an equi- 
librium mass (m + 7n+), an impurity mass m , and the 
system oscillates at the same frequency cujS^ 
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The relation for k' in Eq. (|llaf) is an improvement over 
the low frequency approximation given previously^ Be- 
cause the localization length is analogous to conductivity, 
and the systems described by £ c ^o and £ c ^i occur inde- 
pendently and in parallel, the localization length over all 
values of impurity concentration can be approximated by 
a sum of the two>S 



(12) 



When the concentration of the impurities increases to 
cA > 1, the energy eigenstate becomes more complicated 
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FIG. 3: Localization length £ as a function of impurity con- 
centration c for systems having displacement wavelength A = 
31.8. The three solid curves, from upper to lower, are for a = 
0.30, 0.50, and 0.70. The dashed lines are the dilute limit 
localization length £ c ->o- 



than the dilute impurity example shown in Fig. ^ Ex- 
ample initial conditions for the highest impurity concen- 
tration considered in this experiment (c = 0.5) are shown 
in Fig.[21for three scattering cross sections. At these high 
concentrations, although the frequency remains constant 
everywhere, the wave structure in the disordered region 
is complex. 



G. Parameter Space 

For the nonlinear systems considered, the incident 
wave has wavelength A =31.8, and the anharmonicity 
parameter f3 =1. For all the systems from which energy 
transport is measured, the added mass m+ will have one 
of three values: 6.605, 10.089, and 15.412. With respect 
to a A = 31.8 displacement in a harmonic system, these 
impurity masses have scattering cross sections 0.30, 0.50, 
and 0.70, respectively. 

The initial displacement is a localized mode with lo- 
calization length £. The localization lengths calculated 
from Eq. I|12l) for A = 31.8 and the three impurity masses 
mentioned are shown in Fig. [3] as a function of impurity 
concentration c. The dilute limit result £, c —>o is shown as 
a dashed line, and begins to depart from £ near cA ~ 1. 

Figure 01 also reveals the utility of choosing a nominal 
displacement wavelength A = 32. For cA > 1, the addi- 
tion of impurities has a nonlinear effect on localization 
length. Naturally, one would like to see whether this 
nonlinear behavior has any effect on energy transport. 
For shorter displacement wavelengths, the deviation be- 
tween £ and £ c ^o would not occur until proportionately 
higher impurity concentrations. Alternatively, using a 
longer wavelength would reduce the rate of mode tran- 
sitions considerably, requiring excessively long computa- 
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tional times. 



of energy yields an equivalent alternative expression 
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III. THERMAL CONDUCTION 

There are a number of formal methods for calculat- 
ing thermal conductivity. Thermostats at the bound- 
aries can generate a steady-state flux from which the 
thermal conductivity is calculated using Fourier's law^£ 
The Green-Kubo method 28 is based on time integrals 
of thermal fluxes for a system in thermal equilibrium. 
The Evans NEMD thermal conductivity algorith m 2 ^ 2 ? 
applies a heat field and calculates an averaged heat flux. 
The HelfancU^ moments of the execess energy fluctua- 
tions are calculated for systems in thermal equilibrium. 
In each case, the system is some state of steady-state or 
thermal equilibrium, which does not apply here. 

Elsewhere, the second moment of the site energy has 
been used to characterize pulse dissipationp^iiSiiiiiS 
This approach bears a qualitative resemblance to the 
method of Helfand, which is based on mean squared fluc- 
tuations of excess energy. Here, a brief summary of the 
Helfand moments for thermal conductivity is used to mo- 
tivate the relevance of the second moment of the site en- 
ergy in nonequilibrium chains at zero temperature. 

The thermal conductivity of a collection of freely mov- 
ing particles in thermal equilibrium can be determined 
from energy fluctuations. The energy fluctuation Ei for 
the i-th particle is the difference between the instanta- 
neous site energy Ei and the ensemble averaged value 

(Ei)- 



H e = 
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(13) 



If the energy fluctuation is conserved, and the energy flux 
has a linear dependence on VE, the quantity Ei(x,t) will 
satisfy the diffusion equation. For the boundary condi- 
tions that the initial value Ei(x, 0) is localized about Xio, 
and that £?i(±oo, t) — 0, the solution for Ei(x,t) is Gaus- 
sian, and the measure of spread is the second moment of 



(x — xm) 2 E(x, t) dx ~ 2Ei(x, 0)k£ 



(14) 



where n is the thermal conductivity. 

As Ei represents the fluctuation for a single particle, a 
bulk expression requires an ensemble integral of the sec- 
ond moment. Making no assumption about the indepen- 
dence of particle energies, the thermal conductivity can 
be calculated from a double sum over particle positions^ 



H p = l^Y^ix, - x j{) ) 2 Ei{x,t) Ej{x,Q)j ~2/tf (15) 
Replacing conservation of momentum with conservation 
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Equations (|15|l and l|16f> are the Helfand moments for cal- 
culating the thermal conductivity of a bath of particles. 

There are a number of differences between fluctuations 
in a bath of particles and energy propagation along a 
discrete chain. Equations l|15|) and (|16fl characterize a 
bath of freely moving particles. By contrast, in the FPU 
chain, energy moves, but the masses are, more or less, 
stationary. This is not entirely problematic, however, 
because one can still evaluate the energy that is at Xi. 
The problem can be changed to one in which the energy 
is evaluated at specific points. In this way, the role of 
energy is analogous to concentration in the evaluation of 
self-diffusion. 

Another important distinction is that Eqs. (|15(l and 
(|16[) are functions of the energy fluctuations in a system 
in equilibrium at temperature T. By contrast, the pulse 
moving through the FPU chain is a system that is not in 
equilibrium, and the portion of the chain farthest from 
the initial disturbance is initially at zero temperature. In 
principle, after very long times, the chain would eventu- 
ally reach equilibrium, with the energy distributed over 
all the masses. Because the conceptual problem of inter- 
est is a semi-infinite chain, the equilibrium energy (Ei) 
would approach zero. Under this assumption, the fluc- 
tuation energy E of the bath problem becomes the site 
energy Ei of the non-equilibrium chain problem. 

By analogy to Eqs. (|T5jl and (|TB|l . the energy transport 
in the FPU-/3 chain will be characterized by the second 
moment of the energy. Assuming that the initial pulse 
occupies a small portion of the entire system, a useful 
measure is the second moment about zero: 
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The position — ia is the equilibrium location of the «-th 
mass. The quantity G is an effective transport coefficient 
that is neither self-diffusion nor thermal conductivity. To 
eliminate the effects of fluctuations, the initial value is 
subtracted from the subsequent values. In addition, the 
equation is generalized to allow for arbitrary powers of 
E,: 
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M„(0) ) ~ 2G n t (18) 



These definitions are similar to those used by Frolich et 
al.p^ii and are consistent with the local concept of thermal 
transport of Wagner et ali^ 

A comparison among Mi, M2, H p , and H e is made 
from the early response of a system of length 8000 that 
is initially in a localized mode. The wavelength is 31.8, 
the impurity mass is 11.089, the impurity concentration 



FIG. 4: Wave displacement g; along a chain at various times. 
Chain length is 8000, added mass m+ is 10.089, and impurity 
concentration is 0.010. Each curve represents a time difference 
of 2000, and is offset by a value of one for demonstration 
purposes. Dashed arror denotes ballistic propagation. 
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FIG. 5: Moments Mi, M 2 , H p , H e as a function of time for 
the system shown in Fig. 0] Quantities are normalized, using 
localization length £ and total energy Et, to make the values 
dimensionless. 



FIG. 6: Displacement qi for sinusoid initial condition with 
A = 31.8. 

in Fig. [3] All values are normalized by the initial local- 
ization length (£ ~ 155). H p and H e are also normalized 
by the total energy to be on the same scale as M\ and 
M2. The two pairs of equations correspond well to one 
another. Generally, Mi and M2 are greater in value than 
the corresponding H p and H e . At the shortest times, Mi 
and H p give the nearly the same value. 

H p and H e arc equivalent descriptions of thermal con- 
ductivity. Therefore, the difference between H p and H e 
for the case of an initially localized pulse demonstrates 
that thermal conductivity for these systems is not well- 
defined. As a result, the second moments M n character- 
ize some effective, yet undefined, transport coefficient. 

These curves are also instructive in pointing out the 
distinction between pulse propagation and energy prop- 
agation. From Fig. E| it is clear that low frequency waves 
propagate ballistically through these systems, starting 
from near t — 0. By contrast, Fig. [5] shows no ballistic 
behavior, suggesting that the vast majority of the energy 
is in the higher frequency modes located within the ini- 
tially localized region. Moreover, the moments shown in 
Fig. [5] all continue to increase after there has been suffi- 
cient time for the low frequency ballistic modes to reach 
the far end. 



is 0.010, and /3 = 1. The displacement q along the sys- 
tem is shown in Fig. 01 at time intervals of Aui a t — 2000. 
(The curves are offset vertically from one another, by 
a distance a, for comparison purposes.) The data in the 
figure show that long wavelength displacements move vir- 
tually ballistically along the chain (parallel to dashed ar- 
row) while the higher frequency displacements propagate 
a considerably shorter distance over the same time inter- 
val. 

At regular time intervals, Mi and M 2 are calculated, 
along with H p and H e (assuming that E = E^ and 
(Ei) — 0). The results of the calculations are shown 



IV. INITIAL CONDITION 

The initial condition, composed of a localized mode, 
has the advantage of spatial stability. In the absence 
of anharmonicity, the initial wave will oscillate with a 
periodic amplitude that is constant in time. 

Alternatively, this experiment could have been per- 
formed using either an instantaneous impulse or a short 
sinusoidal pulse like that shown in Fig. |S| Both of these 
initial conditions, however, have drawbacks. Impulses 
will impart relatively little energy to the system unless 
the impulse is large. The sinusoidal wave, not being lo- 
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FIG. 7: Mi as a function of time i for systems with anhar- 
monic parameter /3 = and impurity concentration c = 0.010 
for systems having length L = 8000. 

calized by the impurities, will immediately begin propa- 
gating ballistically along the chain. If the chain had only 
harmonic interactions, the wave would eventually settle 
into a localized state. In a nonlinear chain there will be 
immediate competition between localization of the ini- 
tial wavelength and anharmonic effects that would lead 
to mode transitions that would then experience the same 
competition. 

The time required for the sinusoidal pulse to become 
localized in a harmonic chain is shown in Fig. [7] for the 
A = 31.8 pulse shown in Fig. [fj] The pulse is located 
at one end of the system, the envelope is the hyperbolic 
tangent function, and the initial velocity is zero. The 
harmonic chains have mass impurities located at a site 
with probability c. Three different different values for 
to/ are considered, corresponding to a = 0.30, 0.50, and 
0.70. The data in Fig.[7|suggest that somewhere between 
10 and 100 impurities, depending on the impurity cross 
section, must be encountered before the sinusoidal wave 
becomes localized. The localization length for these three 
cross sections are 298, 155, and 90, respectively, so the 
transient period is approximately ten times the localiza- 
tion length. For the lowest impurity concentrations con- 
sidered in this experiment, this transition length would 
have made the required system length prohibitively long. 



V. EQUIPARTITION 

The rate of both mode transitions and spatial energy 
equipartition will influence the response of the systems. 
The mode transition rate for similar chains configured as 
small hoops initially excited in a single eigenmode stud- 
ied previously^ suggests that impurities initially hasten 
the decay of energy in the excited mode. Over long times, 
however, the rate diminished because the energy became 
localized at the impurities. Because the impurities were 



heavier than the background, the oscillations at the im- 
purities was smaller, reducing the rate of energy loss be- 
cause of the quadratic dependence of amplitude. 

The systems studied here, however, have energy ini- 
tially localized at one end, with the energy already local- 
ized at the impurities. Once transitions start to occur, 
the new modes, which are not localized over the same 
length scale, will propagate and both spontaneously de- 
cay and collide with impurities. 

A. Localization Parameter 

As the mode transitions occur, energy will propagate 
along the chain, redistributing energy. As the energy be- 
comes more evenly distributed among the masses, the 
energy becomes less localized. A measure of how uni- 
formly the energy is distributed among N masses is the 
localization parameter r>i^*i& 

r = n ( ^ E * „ \ (19) 
\(E^) / 

The value of T is a minimum for ergodic behavior and 
increases as the degree of localization increases. 

The localization parameter can be used to estimate the 
number of masses over which energy is distributed. The 
maximum value of T is N, when all the energy is localized 
at one mass. At long time, T approaches a constant, 
Too = T(t — > oo), that only depends on the value of f3m& 
For the FPU-/? system, with = 1, the equilibrium value 
is approximately 1.8. Thus, Too/T is approximately 
equal to the portion of the chain over which energy is 
distributed. 



B. Participating Modes 

Another useful measure of ergodicity is the number 
of harmonic modes contributing to the overall energy at 
a mass. The systems studied will be initially excited in 
one mode (in frequency space). The time required for the 
system to excite the maximum number of modes should 
correspond to the time required for the system to become 
ergodic. 

The energy in a particular mode E w is estimated from 
the harmonic approximation involving the Fourier trans- 
formed (FT) displacement Qi(u>) and momentum Pi{ui) 
of mass m j : 

E ijU = ~ (m^Qj + (20) 

The energy is then normalized using the total number of 
frequency modes N u considered in the FT: 

ei, w =E itU /J2E itU (21) 
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These normalized energies are a measure of energy en- 
tropy Si at mass m^ 15 ! 18 ! 33 ! 34 



Si 



ln(e ijW ) 



(22) 



If all the energy is in a single frequency mode, S equals 
0. If the energy is distributed evenly among all frequency 
modes, S — InN^. 

The alternative expression for the energy entropy is 
exp (S). This quantity is the equivalent number of modes 
contributing to the overall entropy if the energy is uni- 
formly distributed among those modes. To make compar- 
isons among results using different values for N u , results 
are expressed as the fraction of modes n^: 



iW = 



exp {Si 



(23) 



Periodically, nj is calculated at various masses and the 
reported value, n w (f), is the average of these values. 



C. Hoop Example 

The localization parameter T and the fraction of par- 
ticipating modes n w were developed to study systems in 
which the energy is initially distributed throughout. For 
the systems studied here, the initial energy is intention- 
ally localized at one end of the system. If such a system 
was divided into two equal halves, the values for F and 
nu in one half would be very different from the ones cal- 
culated for the other half of the system. Nonetheless, the 
parameters do have utility for these systems. 

One application is the study of behavior within a 
short section of chain. If the section of chain is quasi- 
localized (very few modes are present, oscillating with 
nearly constant amplitude), mode transitions are the pri- 
mary mechanism for inducing energy transport. A hoop 
(periodic boundary conditions), initially excited in one 
wave number mode, could be used to study the behavior 
of a similar section within a much longer section that is it- 
self quasi-localized. The time-dependent behavior of the 
localization parameter T and the relative number of par- 
ticipating modes n w would characterize energy redistri- 
bution, with respect to both space and mode frequency. 

A brief numerical calculation of T and n u is made to 
study the effect. The initial condition is a BID hoop 
having length 318, wavelength 31.8, and (3=1. This 
initial condition differs from a localized state. In a dis- 
ordered system, a single k mode will excite multiple u> 
modes, accelerating the initial rate of mode transitions. 
Nevertheless, the results illuminate general behavior for 
nonlinear BID systems. 

The effects of impurity concentration on energy dis- 
tribution and on mode production are shown in Figs. |HI 
and OH For dilute systems (cA < 1), the time required for 
energy to be distributed among all the masses is nearly a 
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FIG. 8: Localization parameter F as a function of time t for 
a periodic system with length 636, initial wavelength 31.8, 
impurity cross section 0.5, and anharmonicity 1.0. 
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FIG. 9: Fraction of participating modes 71(jj ELS cL function of 
time t for a periodic system with length 636, initial wavelength 
31.8, impurity cross section 0.5, and anharmonicity 1.0. 



constant. For concentrated systems (cA > 1), the spatial 
redistribution of energy is slowed dramatically. 

This behavior is consistent with the data for in 
Fig. The rate that new modes are produced in dilute 
systems is a constant, until n w approaches its asymptotic 
value. At concentrations for which cA > 1, however, the 
initial rate of mode production does not reach the com- 
mon initial rate. The asymptotic value for n u decreases 
with increasing impurity concentration. Therefore, con- 
centrated systems produce fewer modes at a slower rate 
than less concentrated systems. 

The decreasing asymptotic value for !i u with increas- 
ing impurity concentration indicates that modes are sup- 
pressed significantly in concentrated systems. This is 
consistent with the effect of impurities on the spec- 
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tral density u of BID systems. The presence of im- 
purities forces zeroes in u{u>)^^ thereby constrain- 
ing mode transitions. Moreover, as impurity concen- 
tration increases, the spectral density in the interval 
\4K/mj < uj 2 < AK/m ~\ (assuming mj > m D ) becomes 
increasingly suppressed, eventually containing isolated 
delta functionsi2L2& If the frequency of the initial dis- 
placement is in the interval [0 < uj 2 < AK/mi\ , the mode 
will be in a continuous portion of u(oj) for all values of 
impurity concentration. As the impurity concentration 
increases, new modes will be generated more slowly, and 
in proximity to the originial frequency. 

VI. LENGTH AND TIME SCALES 

At the lowest impurity concentrations, the chain is 
composed of long segments of homogeneous nonlinear 
chain. As new modes are generated, these modes are 
not localized and begin to propagate ballistically along 
the chain. The wave will continue to propagate until it 
scatters from an impurity or undergoes a spontaneous 
transition. A spontaneous transition along the homo- 
geneous portion of the chain is unlikely to occur in the 
time required to span the distance between impurities. It 
is more likely that the impurities will initiate scattering 
and mode transitions. The mean free path A = (ccr) -1 
characterizes the effective transport coefficient D(w): 

D(lo) = v g A (24) 

Based on this assumption, thermal dissipation coefficient 
should be proportional to c _1 at the lowest concentra- 
tions. 

As the impurity concentration increases, the localiza- 
tion length decreases and the wave experiences consider- 
ably more wave interference, and transport is more dif- 
fusive because any propagation is now limited by local- 
ization. The relevant time scale is the time required 
for energy to diffuse a distance comparable to the local- 
ization length £: 

DM = L (25) 

k 

When this occurs, the thermal transport coefficient 
should become proportional to c~ 2 . 

VII. RESULTS 

The primary experimental parameter was impurity 
concentration c. The impurity concentration lower limit 
was constrained by computing resources. A previous 
study of these systems revealed that approximately 32 
impurities are required in a chain to ensure reliable en- 
semble statistics^ The upper concentration limit was 
0.5. Above this concentration, the behavior is indistigu- 
ishable from a complimentary study of the pure system 
having mass m; and impurity mass m . 




FIG. 10: The ratio Mi/£ as a function of time for the systems 
having impurity concentration 0.500 and length 16 000. The 
j3 = 1 data have positive slopes and all the f3 = data fall on 
top of one another. 



Some of the results include estimates of uncertainty for 
quantities calculated from ensemble averages. For a cal- 
culation performed on an ensemble of W systems, there is 
a population standard deviation s and an average value. 
For this study, the average value is the meaningful quan- 
tity. The uncertainty in the reported average value is 
s I y/W, and is referred to here as the standard deviation 
in the mean (SDM). 

There are two special cases within the parameter space 
that require special consideration. The M n data for the 
nonlinear systems are only meaningful when the M n data 
for the corresponding harmonic system is a constant. 
This must be true for all cases, especially at concentra- 
tions for which cA 3> 1. 

In theory, M n for harmonic systems would be a con- 
stant for all time. In practice, the mapping of the con- 
tinuum system onto the discrete lattice, and the reloca- 
tion of the ends, introduced a small amount of instability 
that led to small fluctuations in M n . These fluctuations, 
however, were far smaller than the changes for the an- 
harmonic systems. 

As a brief example, Fig.^|is a plot of Mi as a function 
of time for systems having having impurity concentration 
0.500 and length 16 000. In the figure, the = 1 data ap- 
pear as lines having positive slope. The harmonic = 
data for all three impurity masses lie upon one another 
near M\ = 0. Figure is doubly instructive. It demon- 
strates that the oscillations for the = data are negli- 
gible, even for the systems for which Mi has the smallest 
values. Moreover, the harmonic data remain stable in 
concentrated systems: cA 3> 1. 
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TABLE I: Interval from which S n was calculated for systems 
having length L. 



L 


Interval 




16000 


5000 < loo t < 


20000 


32000 


10000 < L0 o t < 


25000 


64000 


20000 < loo t < 


50000 


96000 


20000 < L0 o t < 


80000 
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FIG. 11: The time exponent <5i for Gi as a function of im- 
purity concentration c for impurity cross sections a = 0.30, 
0.50, 0.70. Numbers between vertical dashed lines denote N 
(k=1000). 
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FIG. 12: The time exponent 82 for G2 as a function of impu- 
rity concentration c for impurity cross sections a — 0.30, 0.50, 
0.70. Numbers between vertical dashed lines denote system 
size used (k=1000). 



than for Si. Moreover, it was difficult to establish a pre- 
cise value for 62 at higher concentrations. As a result, 
subsequent analysis is confined to Mi. 

According to Li et al.^ if the mean squared displace- 
ment of a particle is proportional to t a with (0 < a < 2), 
the thermal conductivity can be expressed as a function 
of the system size L: k ~ L b with 6 = 2 — 2 /a. There- 
fore, a system obeying normal diffusion implies normal 
heat conduction obeying Fourier's law (6 = 0). 



A. Time Exponent 

A quantitative characterization of M n is made by as- 
suming a power-law dependence on time t: 

M n = 2G n t 5 " (26) 

The first task is to determine the value of the exponent 
so as to distinguish among ballistic, diffusive, or sub- 
diffusive transport. The calculation of 5 n is affected by 
initial transients, and the details of the analysis are given 
in the Appendix. The time interval from which the val- 
ues of S n were calculated were constrained by the initial 
transients at small times and total energy conservation 
at long times. The specific intervals are shown in Tabled 
for each system length. 

The results of the analyses for <5i are summarized in 
Fig.^Jfor all the systems. (In the figure, the symbols at 
a particular concentration are displaced horizontally to 
distinguish individual error bars representing the SDM.) 
With only one exception, the estimated values for 5\ are 
within one SDM of 1. Therefore, subsequent analysis of 
M\ is based on linear model with respect to time. 

Results of the analysis for 82 are shown in Fig. 1121 for 
the same systems. Although the values of 62 are near 1 
for most systems, there is considerably more variability 



B. Transport Coefficient 

Because Si was shown to be within one SDM of 1, the 
determination of 2Gi is based on the assumption of a 
linear relationship between Mi and t over the same in- 
tervals shown in Table [I] Although Mi was determined 
by linear regression with a linear model, calculating the 
uncertainty in Mi required a slightly more involved anal- 
ysis; details are give in the Appendix. 

Estimates of 2Gi for all the systems considered are 
plotted in Fig. ^| as a function of the impurity concen- 
tration c. The calculated values appear as filled symbols 
having error bars that represent the SDM. The thin solid 
lines connecting the symbols are only to guide the eye. 
The two dashed line segments appearing above the data 
indicate slopes that are proportional to 1/c and 1/c 2 . 
The three dotted lines are proportional to £, as they ap- 
pear in Fig. El 

From the data in Fig. 1131 there appear to be three 
regions of interest. At low concentration, the transport 
coefficient Gi is proportional to c , as was expected 
from the discussion of time and length scales. As the 
impurity concentration increased, and the diffusion time 

decreases sufficiently that Eq. lj25H dominates, and Gi 
becomes proportional to c~ 2 . This transition is not ap- 
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FIG. 13: Transport coefficient G as a function of impurity 
concentration c for impurity cross sections a = 0.30, 0.50, 
0.70. Dashed line is proportional to c _1 . Dotted line is pro- 
portional to £. 



parent in the results of Payton et alw^ on a similar system 
having thermostats because their computing resourses 
prevented them from resolving the smaller concentrations 
required to see the effect. 

There is a noticeable difference in the rate at which 
the three curves make the transition from c _1 to c~ 2 . 
Based on the transition data for similar FPU-/3 systems, 
the transition rate should be inversely proportional to 
the impurity mass m/^i Therefore, the transition rate 
for the m + = 15.412 is slowing more rapidly, for a given 
change in impurity concentration, than the other two im- 
puriy masses. As a result, the m+ = 15.412 systems do 
not exhibit distinct 1 [c dependence at the concentrations 
studied. 

The interesting behavior was the final transition at 
higher concentrations to a transport coefficient that is 
proportional to the initial localization length £. The 
three dotted curves labelled oc £ are proportional to the 
localization lengths that appear in Fig. [3] for A = 31.8. 
The value of £ is multiplied by the same coefficient (ap- 
proximately 0.09) to make the curves best agree with the 
measured transport coefficient. 



C. Localization Parameter 

The localization parameter can also be used as a mea- 
sure of energy propagation. Although systems of differ- 
ent lengths were used, energy transport, starting from 
a localized state, should be independent of total system 
length. Based on the previous discussion of the localiza- 
tion parameter, the ratio L/T is proportional to the num- 
ber of masses over which the total energy is distributed. 
For a given mi and c, and assuming that L>(, the num- 
ber of masses, and not the fraction of the system (Too/T), 
should also be independent of total system length. Thcrc- 
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FIG. 14: Localization parameter Y as a function of time t 1 ' 2 
for systems having m+ = 6.605. 
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FIG. 15: Localization parameter T as a function of time t 1 ^ 2 
for systems having m + — 10.089. 



fore, the ratio L/T will serve as a means for comparing 
results from systems of different lengths. 

If the linear time dependence of the second moment 
M\ indicates diffusive behavior, the number of masses in- 
volved in energy transport should increase in proportion 
to t 1 / 2 . In addition, because transport in dilute systems 
is proportional to c _1 , the quantity cLr should approach 
a constant for systems with low impurity concentrations. 
The quantity cL/T, as a function of i 1 / 2 , is plotted in 
Figs. H4H16l for most of the systems studied. In the fig- 
ures, curves for smaller impurity concentration appear 
consecutively lower in the graph. 

There are two noteworthy features in Figs. ll4Hlb1 The 
number of masses participating appears to be linear over 
the time intervals considered, particularly for the dilute 
impurity systems. This is consistent with the expecta- 
tion of diffusive energy transport. Also, as expected, the 
curves for the most dilute impurity concentrations appear 
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FIG. 16: Localization parameter F as a function of time t 1 ^ 2 
for systems having m+ = 15.412. 
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FIG. 17: Transport coefficient G as a function of /3 for two 
impurity concentrations: c = 0.100, 0.005, and a — 0.70. 



to approach an asymptote, supporting the c _1 depen- 
dence for the number of masses participating in energy 
transport. 



VIII. DISCUSSION 
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FIG. 18: The ratio M\(t m ax)/£ 2 as a function of impurity 
concentration c. 



of (3. Over the time interval studied, values of G\ de- 
creased for values of (3 below 0.5, falling to near zero at 
(3 = 0.2. For 0.5 < (3 < 5.0, the transport coefficient has 
a weak dependence on the anharmonicity, which is con- 
sistent with numerical experiments on finite temperature 
nonlinear BID systems^ Therefore, while the calculated 
values of G are not constant, there does not appear to be 
any significance to a particular value for (3 that is greater 
than 0.5. 



B. Diffusion Time t% 

The coefficient 2G\ in Fig. ^| represents the trans- 
port coefficient that characterizes the time defined in 
Eq. 



(27) 



At low concentration, because both G\ and £ are pro- 
portional to c _1 , the diffusion time must also be pro- 
portional to c _1 . At intermediate concentration, G\ is 
proportional to c -2 , so is a weak function of impurity 
concentration. At the highest concentrations studied, G\ 
and are both proportional to £. At these concentra- 
tions, however, £ is a weak function of impurity concen- 
tration, so the same must be true for t^. 



A. /9-Dependence 

It was argued that the results given for (3=1 are in- 
icative of results for 'large' values of (3 that are above 
the critical threshold that leads to ergodic behavior. As 
a check, the value of G\ for two systems are calcu- 
lated for different values of (3. The two systems are 
(c = 0.005, cr = 0.70) and (c = 0.100, cr = 0.70). The 
analysis for these systems was carried out in a manner 
identical to that for the (3—1 data. 

The values for G\ are shown in Fig. 1171 as a function 



C. Time Scale 

The significance of the reported results presented de- 
pends, in part, on the relative duration of the calcula- 
tion. One measure of duration, to gauge whether long 
times have been probed, is the ratio of the distance en- 
ergy propagates along the system to its initial span. The 
maximum energy propagation length is proportional to 
Mi(t max ), and the initial span is proportional to the ini- 
tial localization length £. 
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The ratio Mi(t max ) / £ 2 for each system is plotted in 
Fig- El The ratio varied from 10 to 100, increasing with 
impurity concentration. The square root of this ratio is a 
measure of the depth to which energy propagated to the 
initial depth energy was distributed. Therefore, energy 
penetration depth, relative the initial localization length, 
varied from 3 to 10. 

Extending this experiment to longer times will prob- 
ably yield similar results. It is unlikely that ballistic 
transport will occur at long times. As the calculation 
progresses, the energy is being distributed over a greater 
number of masses. Because both energy and transition 
rates are proportional to the square of the oscillation am- 
plitude, the mode transition rate is deceasing in propor- 
tion to t^ 1 / 2 



D. Memory Effect 

Although previous arguments address why the trans- 
port coefficient should vary from c _1 to c -2 , they do not 
explain why the transport coefficient should be propor- 
tional to £ at the highest concentrations. More specifi- 
cally, it was unexpected that the transport coefficient is 
proportional to the original localization length and not 
some averaged value. Figure EJ shows that the energy 
has propagated nearly ten times the localization length 
for the highest impurity concentrations, yet these systems 
seem to retains some "memory" of the initial state. 

An explanation may be made from the results of the 
hoop experiment shown in Fig. For those systems, 
both the number of modes and the transition rate de- 
creased significantly with increasing impurity concentra- 
tion. At these concentrations and time scales, mode gen- 
eration may have been dominated by impurity scatter- 
ing. As the localization length is quite small for most 
modes, any new mode could only propagate a short dis- 
tance before becoming localized. If the impurities con- 
strained mode generation sufficiently so that generated 
modes had frequencies near the initial mode, the initial 
localization length would remain as the controlling length 
scale over which new modes could propagate before scat- 
tering again. 



IX. CONCLUSION 

For systems initially localized near one end of a FPU- 
(3 chain, anharmonicity facilitates energy transport along 
the chain. For energy eigenstates of harmonic disordered 
chains, sufficient anharmonicity leads to diffusive energy 
transport. The second moment of the site energies was 
linear over times long enough for energy to have propa- 
gated a distance nearly equal to ten times the initial spa- 
tial extent of the pulse. Although direct mass displace- 
ment observation showed low frequency modes propagat- 
ing ballisticall, the energy content in these modes was 
insufficient for super-diffusive energy transport. 



Alternate measures of energy and frequency distribu- 
tion corroborated the assertion of diffusive behavior. Wc 
found that the localization parameter was proportional 
to the square root of time, suggesting that energy dis- 
tributes itself among masses in a diffusive manner. This 
diffusive short range behavior is consistent with the ob- 
servation of more long range diffusive behavior. Cal- 
culations of the relative number of modes participating 
in periodic (hoop) systems suggest that an increase in 
the number of impurities reduces the number of active 
modes. 

The most interesting aspect of the transport coefficient 
was the concentration dependence. At the lowest im- 
purity concentration c values, the transport coefficient 
was proportional to c . At higher impurity concentra- 
tions, the transport coefficient developed a c~ 2 depen- 
dence. These two concentration dependences are consis- 
tent with arguments based on the competition between 
diffusion times and mode transition rates. 

At the highest impurity concentrations, the transport 
coefficient was proportional to the original localization 
length. The postulated cause for this was from the se- 
vere suppression of new frequency modes that occurs in 
concentrated BID systems. These new modes, having a 
frequency in proximity of the original frequency, experi- 
enced a similar localization length, which dominated the 
distance travelled before subsequent scattering events. 
As a result, the original localization length remained the 
dominant length scale over which transport occured for 
any given phonon. 
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APPENDIX: DATA ANALYSIS 

The analyses of the results involve some minor sub- 
tleties that deserve a clear exposition. Neglecting these 
subtleties and performing an ordinary least squares 
(OLS) analysis of the data would lead to misleading re- 
sults. Specifically, not adjusting for initial transient be- 
havior would lead to a different conclusion regarding the 
existence of diffusive energy transport. 

As mentioned previously, the reported uncertainties 
are the standard deviation in the mean (SDM) calculated 
from the ensemble population standard deviation s. For 
an ensemble of W systems, the SDM reported here is 
s/y/W. This uncertainty characterizes the reported av- 
erage value from the population. 
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FIG. 19: Mi /£ 2 as a function of time. System length is 64 000, 
impurity concentration is 0.002, single impurity cross section 
is 0.30. The measured data appear as a solid line, the shifted 
data appear as dashed line, and the line segment denotes the 
range over which regression was performed. 

1. Energy Fluctuation 

For each calculation, the total energy T(t) was cal- 
culated as a function of time. Due to randomness, the 
initial total energy fluctuates among the ensembles. To 
ensure that values of T were on comparable scales, the 
values were divided by the initial value T(0). The ratio 
Y m (i)/T m (0) (1 < m < W) i is calculated as a function 
of time for each of the W systems composing the ensem- 
ble. The averaged values are calculated at each of the P 
values of U : 

m) = JL y i<i<p (a.i) 

m v ' 

These averaged values, along with the population stan- 
dard deviation, are pooled and stored in the output data 
file. The population of W values at fj do not, unfortu- 
nately, represent energy fluctuation for a single system. 
Rather, it characterizes statistical fluctations among the 
different systems, evaluated at the same time. 

Energy fluctuation can only be approximated from flu- 
cuations in the P average values of T(ti). The standard 
deviation Sy of each T(ti) represents a standard devia- 
tion in the mean. The population standard deviation is 
approximated by s^VP- This population standard de- 
viation was never more than 0.2 %. Spot checks of Y(t) 
in individual systems rarely gave a standard deviation 
greater than 0.2 %. 

2. 5 n Analysis 

The straightforward means of determining the time ex- 
ponent S n is from the slope of a log-log plot of M n versus 
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FIG. 20: A log-log plot of data in Fig. ED along with hori- 
zontal error bars denoting the population standard deviation; 
the vertical risers are omitted for clarity. The solid line was 
determined by OLS regression applied to the mean values. 



time t. This approach assumes that the power law de- 
pendence exhibits itself at t — 0. In reality, there is a 
transient period, after which the value of M n begins to 
increase. Although this increase appears to be linear, 
the data must be corrected to eliminate the effect of the 
transient. 

The effect of the transient behavior is negated by shift- 
ing the data for M n vertically. OLS linear regression 
is applied to the M n versus t data over the time inter- 
vals given in Table |U The intercept calculated from the 
OLS regression is subtracted from M n (t). The process is 
demonstrated in Fig. ^| for one particular system. The 
original Mi data are shown as a solid line. Regression is 
applied to the interval [20000,50000], and the data are 
shifted vertically so that the linear approximation to this 
interval has zero intercept. 

OLS linear regression is then applied to the logarithm 
of these adjusted data versus the logarithm of time. The 
slope of this regression calculation is the value reported 
in Figs. Illlandll2lfor Si and 62, respectively. 

The uncertainty in 5 n is calculated from both the re- 
gression residuals and the ensemble population of M n 
values. Fortunately, because the population standard de- 
viation in M n increases in time, the logarithm transform 
yields uncertainties that are nearly constant over the time 
intervals of interest. The adjusted data from Fig.lTTflare 
shown in Fig. I20I along with the ensemble population 
standard deviations that are denoted by horizizontal er- 
ror bars; the vertical risers are omitted for clarity. (The 
ensemble population standard deviation, instead of the 
SDM, is shown in the figure for clarity of the demon- 
stration. For an ensemble of W systems, the SDM error 
bars would be a factor of \fW smaller than those shown 
in the figure.) Also shown in the figure is the result of 
the regression analysis of the average values. The error 
bars denote the ensemble uncertainty, and the residuals 
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represent the regression uncertainty. 

The OLS linear regression in log-space will yield an 
estimated uncertainty (standard deviation) for slope that 
is a function of the standard error of the residuals. The 
standard error s can be used to estimate the regression 
standard deviation s reg for the slope S n - 

(A.2) 



V &XX 

where the quantity Sxx is the sum of squares: 

, 2 



(A.3) 



The quantity hit represents the average value over the 
specified interval. 

The uncertainty in S n should also reflect the SDM 
recorded for the ensemble s ens . It is assumed that these 
two uncertainties are independent of one another, and 
that they are additive. The uncertainty (estimated stan- 
dard deviation) in the time exponent is 



2 ^res 

S S, 



S 2 s 2 



Sxx Sxx 



(A.4) 



Because s ens is the majority of sg n , the uncertainty in 
S n is referred to as a SDM, and the coverage factor is 
approximately equivalent to one standard deviation for a 
normal distribution. 



3. Mi Analysis 



To use OLS techniques, the error term must be additive 
and constant. Equation (|A.5|I can be transformed into a 
suitable model: 



M n (U) _ A 



2G n + a 



(A.6) 



Unfortunately, the initial transient behavior that con- 
founded the calculation of S n must also be accounted 
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FIG. 21: Transformed Mi to demonstrate how uncertainty in 
Si is calculated. 



Based on the results of the Si analysis, a hypothesis 
of a linear relationship between the second moment Mi 
and time t cannot be rejected. The ability to use a lin- 
ear model simplifies the analysis of M\. The transport 
coefficient 2G\ is the slope that is calculated from OLS 
regression using a linear model. 

The uncertainty in 2G\ cannot be determined from a 
regression analysis of the Mi versus t data because the 
uncertainties in Mi increase in time. Fortunately, the 
data can be transformed into a more suitable format. 

Assuming diffusive behavior, the values of Mi are a 
collection of lines, radiating from the origin. The error 
model for the observations assumes that there exists an 
inherent error e and that the total error increases with 
time: 

M n (ti) = A + 2G n t i + t i e i (A.5) 



for here. The error model assumes that the uncertainty 
grows linearly from t = 0. In reality, the increase in M n 
occured after some initial transient time t . To correct 
for this, linear regression is applied to the Mi versus t 
data to determine the value of t . Using to, the data 
are shifted horizontally so that the linear region of in- 
terest points to the origin. These shifted data are then 
transformed according to Eq. I|A.6|) . 

The result of this transform, applied to both the mean 
and SDM, is shown in Fig. [21 for the data in Figs. IT§1 
and [23 In Fig. [21 the error bars represent the SDM. 
Over the range of regression (see Table HJ, the SDM is 
relatively constant, consistent with Eq. 1|A.6|) . Because 
the model in Eqn. (|A.6(1 assumes a constant factor of 2G\ , 
the error bars in Fig. [21 represent the SDM for 2G\ . 
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